
// Created by me1on on 22-5-17.

#include "RefinedRectangleReceiverRectangleGridRayTracing.cuh"
#include <algorithm/rayTracing/dataLoader/CloudLoader.h>
#include <algorithm/rayTracing/rayTracingDataStructure/receiver/CudaRefinedRectangleReceiver.cuh>
#include <algorithm/rayTracing/rayTracingPipeline/quasiMonteCarloRayTracing/Grid/rectGridDDA.cuh>
#include <algorithm/shared/tool/randomNumberGenerator/RandomGenerator.cuh>

using namespace solar;
using namespace solar::cuda;

__global__ void map_raytracing(SunrayArgument sunrayArgument,
                               CudaRefinedRectangleReceiver refinedRectangleReceiver,
                               CudaRectangleGrid rectangleGrid, HeliostatArgument heliostatArgument,
                               float3* d_subheliostat_vertexes, float factor, bool res,
                               CloudLoader cloud)
{
    int myId = global_func::getThreadID(); // 线程ID
    int numberOfMicroHeliostats =
        heliostatArgument.micronums * heliostatArgument.helionums; // 所有微定日镜个数
    if (myId >= numberOfMicroHeliostats * sunrayArgument.numberOfLightsPerGroup)
    {
        return; // id比光线数目大，返回
    }
    int microsize = heliostatArgument.micronums; // 每个定日镜的微定日镜个数
    int no = myId / microsize / sunrayArgument.numberOfLightsPerGroup; // 第几号镜子

    int helio_no =
        myId / sunrayArgument.numberOfLightsPerGroup / heliostatArgument.micronums; // 第几个定日镜
    int micro_no = myId / sunrayArgument.numberOfLightsPerGroup %
                   heliostatArgument.micronums; // 该定日镜的第几个微定日镜

    int substart = helio_no * 3;

    float3 helio_pos =
        (d_subheliostat_vertexes[substart] + d_subheliostat_vertexes[substart + 2]) / 2;
    float3 helio_norm = normalize(
        cross(d_subheliostat_vertexes[substart] - d_subheliostat_vertexes[substart + 1],
              d_subheliostat_vertexes[substart + 2] - d_subheliostat_vertexes[substart + 1]));

    int col_no = micro_no / heliostatArgument.col;
    int row_no = micro_no % heliostatArgument.col;
    float3 micropos = helio_pos;
    float3 micronorm = helio_norm;

    micropos.z = col_no * heliostatArgument.pixel + heliostatArgument.pixel / 2 -
                 heliostatArgument.sizex / 2;
    micropos.y = heliostatArgument.sizey / 2;
    micropos.x = row_no * heliostatArgument.pixel + heliostatArgument.pixel / 2 -
                 heliostatArgument.sizez / 2;

    float3 local = micropos;
    local = global_func::local2world_rotate(local, micronorm);
    local = global_func::translate(local, helio_pos);
    micropos = local;

    int address = (myId) % sunrayArgument.pool_size;

    // 采样数组的采样值
    float3 dir = global_func::local2world_rotate(sunrayArgument.d_samplelights[address],
                                                 -sunrayArgument.sunray_direction);
    // 入射光方向（世界坐标系）
    // dir = -sunrayArgument.sunray_direction;
    // 计算第helio_no个定日镜的第micro_no个微定日镜的法向和坐标

    float3 origin = micropos; // 微定日镜中心坐标
    dir = normalize(dir);
    if (rectGridDDA::collision_(origin, dir, rectangleGrid, d_subheliostat_vertexes,
                                heliostatArgument, no))
    {
        return;
    }
    float3 normal = micronorm; // 微定日镜法向方向，局部坐标系
    address = (myId) % sunrayArgument.pool_size;
    float3 normal1 = global_func::local2world_rotate(sunrayArgument.d_perturbations[address],
                                                     normal); // 微定日镜法向方向，世界坐标系
    normal1 = normalize(normal1);

    //     //先和其他定日镜求交,看是否被遮挡
    //     //如果没有被遮挡,求光线与云的交点,如果在数组中,取数组中的值,不在,取1.返回值是float值
    //     float3 dir_n= make_float3(0.0,1.0,0.0);
    //     float cos_theta= dot(dir,dir_n);
    //     float len=cloud.start_point.y/cos_theta;
    //     float3 dirs=dir*len;
    //     float2 M= make_float2(origin.x+dirs.x ,origin.z+dirs.z );
    //     float times=cloud.start_point.y/-sunrayArgument.sunray_direction.y;
    //     float2 cloudcenter=
    //     make_float2(-sunrayArgument.sunray_direction.x*times,-sunrayArgument.sunray_direction.z*times);
    //     float2
    //     cloudstart=make_float2(cloudcenter.x+cloud.start_point.x,cloudcenter.y+cloud.start_point.z);
    //     float my=M.x-cloudstart.x;//M.x, M.y 是M点在世界坐标系的（x,z)坐标值
    //     float mx=M.y-cloudstart.y;//需要转换成 cloud map局部坐标
    //     /*世界坐标  --> X
    //      *         |
    //      *         |
    //      *         Z
    //      * cloud map 局部坐标  --> y
    //      *                    |
    //      *                    |
    //      *                    x
    //      * M.x-cloud.start_point.x M点在局部坐标系的第几列
    //      * M.y-cloud.start_point.z M点在局部坐标系的第几行
    //      * A    B
    //      *    M
    //      * D    C
    //      * 根据A,B,C,D四点的值双线性插值得到M的值
    //      */
    //     int Ax= floor(mx/cloud.pixel);//A在cloud map中第几行
    //     int Ay= floor(my/cloud.pixel);//A在cloud map的第几列
    //     float decay=1.0;
    //     if(Ax<=0||Ax>cloud.x_size||Ay<0||Ay>cloud.y_size){
    //         decay=1.0;
    // //      printf("times %f cloudstart %f %f %f %f \n",times,cloudstart.x,cloudstart.y,mx,my);
    //     }
    //     else{
    // //        printf("my %f mx %f\n",my,mx);
    //         /*世界坐标  --> X
    //          *         |
    //          *         |
    //          *         Z
    //          * cloud map 局部坐标  --> y
    //          *                    |
    //          *                    |
    //          *                    x
    //          * M.x-cloud.start_point.x M点在局部坐标系的第几列
    //          * M.y-cloud.start_point.z M点在局部坐标系的第几行
    //          * A    B
    //          *    M
    //          * D
    //          * 根据A,B,C,D四点的值双线性插值得到M的值
    //          */

    //         int A_pos=Ay+Ax*cloud.x_size;

    //         float Q11=cloud.d_cloudMap[A_pos];                  //A
    //         float Q12=cloud.d_cloudMap[A_pos+1];                //B
    //         float Q21=cloud.d_cloudMap[A_pos+cloud.x_size];     //D
    //         float Q22=cloud.d_cloudMap[A_pos+cloud.x_size+1];   //C
    //         //A x1,y1
    //         //B x1,y2
    //         //C x2,y2
    //         //D x2,y1
    //         float x1=Ax*cloud.pixel;
    //         float x2=(Ax+1)*cloud.pixel;
    //         float y1=Ay*cloud.pixel;
    //         float y2=(Ay+1)*cloud.pixel;

    //         decay=(Q11*(x2-mx)*(y2-my)+Q21*(mx-x1)*(y2-my)+Q12*(x2-mx)*(my-y1)+Q22*(mx-x1)*(my-y1))/cloud.pixel/cloud.pixel;
    //         decay=1.0-decay;
    //     }

    dir = normalize(reflect(-dir, normal));

    if (rectGridDDA::collision_(origin, dir, rectangleGrid, d_subheliostat_vertexes,
                                heliostatArgument, no))
    {
        return;
    }

    // Step 3: intersect with receiver and calculate energy.
    float t;
    float u;
    float v;
    if (!refinedRectangleReceiver.GIntersect(origin, dir, t, u, v))
    {
        return;
    }
    float eta = 0.0;
    if (t <= 1000.0F)
    {
        eta = 0.99331F - 0.0001176F * t + 1.97F * (1e-8F) * t * t;
    }
    else
    {
        eta = std::exp(-0.0001106F * t);
    }
    // float energy=fabsf(dot(dir, normal1)) * eta* factor*decay;
    float decay = 1.0;
    float energy = std::fabs(dot(dir, normal1)) * eta * factor * decay;

    // Step 3: Add the energy to the intersect position
    // Intersect position
    int2 row_col = make_int2(u * refinedRectangleReceiver.getResolution().y,
                             v * refinedRectangleReceiver.getResolution().x);
    address = row_col.x * refinedRectangleReceiver.getResolution().x +
              row_col.y; // col_row.y + col_row.x * resolution.y

    //    printf("row_col %d %d\n",row_col.x,row_col.y);
    float* image = refinedRectangleReceiver.getDeviceImage();
    atomicAdd(&(image[address]), energy); // CUDA atomic
}

void solar::cuda::RefinedRectangleReceiverRectangleGridRayTracing(
    SunrayArgument* sunrayArgument, CudaRefinedRectangleReceiver* refinedRectangleReceiver,
    CudaRectangleGrid* rectGrid, HeliostatArgument* heliostatArgument,
    float3* d_subHeliostat_vertexes, float factor, bool res, CloudLoader* cloud)
{
    int n_threads = 512;
    dim3 n_blocks;
    global_func::setThreadBlocks(n_blocks, n_threads,
                                 heliostatArgument->helionums * heliostatArgument->micronums *
                                     sunrayArgument->numberOfLightsPerGroup,
                                 true);

    printf("num of rays= %d \n", heliostatArgument->helionums * heliostatArgument->micronums *
                                     sunrayArgument->numberOfLightsPerGroup);
    printf("nBlocks = %d %d %d, nThreads = %d\n", n_blocks.x, n_blocks.y, n_blocks.z, n_threads);

    //    printf("%lld\n", heliostatArgument.numberOfMicroHeliostats *
    //    sunrayArgument.numberOfLightsPerGroup);
    map_raytracing<<<n_blocks, n_threads>>>(*sunrayArgument, *refinedRectangleReceiver, *rectGrid,
                                            *heliostatArgument, d_subHeliostat_vertexes, factor,
                                            res, *cloud);

    cudaDeviceSynchronize();
    checkCudaErrors(cudaGetLastError());
}
